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In the context of metric perturbation theory for non-spinning black holes, extreme mass ratio 
binary (EMRB) systems are described by distributionally forced master wave equations. Numerical 
solution of a master wave equation as an initial boundary value problem requires initial data. 
However, because the correct initial data for generic-orbit systems is unknown, specification of trivial 
initial data is a common choice, despite being inconsistent and resulting in a solution which is initially 
discontinuous in time. As is well known, this choice leads to a "burst" of junk radiation which 
eventually propagates off the computational domain. We observe another potential consequence 
of trivial initial data: development of a persistent spurious solution, here referred to as the Jost 
junk solution, which contaminates the physical solution for long times. This work studies the 
influence of both types of junk on metric perturbations, waveforms, and self-force measurements, 
and it demonstrates that smooth modified source terms mollify the Jost solution and reduce junk 
radiation. Our concluding section discusses the applicability of these observations to other numerical 
schemes and techniques used to solve distributionally forced master wave equations. 

PACS numbers: 04. 25. Dm (Numerical Relativity), 02.70.Hm (Spectral Methods), 02.70.Jn (Collocation 
methods); AMS numbers: 65M70 (Spectral, collocation and related methods), 83-08 (Relativity and gravi- 
tational theory, Computational methods), 83C57 (General relativity, Black holes). 



I. INTRODUCTION 

Extreme mass ratio binary (EMRB) systems are typi- 
cally comprised of a small compact object, such as a stel- 
lar black hole, orbiting a super-massive blackhole, and 
the gravitational radiation generated by such systems is 
potentially detectable by the LISA project. A number 
of approaches attempt to model the resulting gravita- 
tional waveforms, including effective one body formula- 




tions [!H3j|, effective field theory techniques 
Newtonian expansions 3 , self- force effects 
different gauge choices [ll| - [l3| . When including high- 
order effects or performing comparisons between tech- 
niques, improved EMRB modeling will increasingly re- 
quire the identification and reduction of all error sources 
(both numerical and systematic). 

Consider a small perturbation h^ v of a fixed back- 
ground Schwarzschild metric, where h^ v satisfies the lin- 
earized Einstein equations. The metric perturbation h^ v 
describing an EMRB can be reconstructed from a collec- 
tion of scalar master functions $ , each of which obeys 
a forced wave equation of the form (with all multipole 
indices suppressed) 



- <9 f 2 * + <9 2 * - V(r)V 

= f(r) [G(t, r)S(r - r p (t)) + F(t, r)5'(r - r p (t))] . 



(1) 
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The coordinates here are the areal radius r, the Regge- 
Wheeler tortoise coordinate x = r + 2Mlog(ir/M — 
1), and the time-dependent radial location r p (t) of the 
smaller mass or "particle" . M is the mass parameter of 
the background solution, f(r) = 1 — 2M/r, and V(r) 
is either the Reggc- Wheeler or Zerilli potential (explicit 
expressions for both are given in Sec. Ill B|) . The dis- 
tributional inhomogeneity on the right-hand side of (fT]) 
involves Dirac delta functions, as well as the ordinary 
functions F(t,r) and G(t,r). For all possible choices 
of the master function, F(t, r) and G(t,r) are listed in, 
for example, in Refs. Here it suffices to note 

that their evaluation requires knowledge of the parti- 
cle's four- velocity u a , orbital energy and angular mo- 
mentum parameters (E p , L p ), and equatorial location 
(r p (t),ir/2,<fi p (t)). In the model we study, integration of 
the geodesic equations determines the timelike particle 
trajectory (r„(t), <j> v (t)) in the equatorial plane 9 = tt/2. 

[3 EMI 

One approach for computing EMRB waveforms is to 
numerically solve Eq. (fTJ) as a time-domain initial value 
problem with prescribed initial data. The exact initial 
data for generic point-particle trajectories is non-trivial, 
and the most common choice is therefore to set both ^ 
and its time derivative to zero. (See Refs. (l9l - [22| for 
the construction of more realistic data.) Inspection of 
(TT]) shows that trivial data is inconsistent with the jump 
conditions stemming from the delta function terms in the 
inhomogeneity. As a result, trivial data results in an im- 
pulsive (i. e. discontinuous in time) start-up. This paper 
addresses the main question of if, and when, a physi- 
cal solution eventually emerges from such trivial initial 



2 



data. Ideally, we would have both the correct source 
terms and initial conditions. Without the exact initial 
data, we consider modifying the source terms such that 
they are consistent with the choice of trivial initial data. 
Precisely, the source terms are "switched on" smoothly 
via the following prescription: 

F(t,r) -> F{t,r)x 

f §[erf(\/tf(i - t - t/2) + 1] for to < t < to + r 
\ 1 for t > t + t, 

(2) 

and the same for G(t,r). Typically, the initial time to — 
0, and the timescale r is much shorter than the final 
time of the run. Choosing suitable r and 5, one achieves 
smooth and consistent start-up to machine precision. 

To appreciate some of the issues associated with the 
main question above, consider a particle in a fixed circu- 
lar orbit. The energy Eqw and angular momentum Lqw 
luminosities for gravitational waves are then constant in 
time and obey the relation Eqw = ttLawi where fl is 
the angular velocity of the particle. However, verification 
of this relationship is limited by a finite computational 
domain, leading to an 0(r~ r ) error (see Ref. [23| for a 
recent suggestion towards overcoming this limitation). 
Therefore, numerical verification of Eqw = QLqw is 
a useful diagnostic only in the distant wave-zone. In 
the near-zone we might also test "-Egw = QLqw" , now 
constructing the luminosities with self-force quantities 
via (|26|) below; however, because \£' is discontinuous at 
the particle location, self-force measurements will involve 
large errors unless due care is taken. For generic quasi- 
periodic orbits, selection of a meaningful set of diagnos- 
tics is not straightforward. In particular, we can nei- 
ther infer steady-state behavior throughout the computa- 
tional domain, nor claim we have a solution which solves 
the hypothetical "true" initial value boundary problem. 
These difficulties are due to the inconsistent initial con- 
ditions. That is, we are really solving a problem dif- 
ferent from the physical one. As a partial resolution of 
these issues, we examine a direct test condition which is 
necessary to claim that a physically correct solution has 
been achieved everywhere in the computational domain. 
This is a simple self-consistency condition relating the 
Cunningham-Pricc-Moncrief (CPM) and Rcggc- Wheeler 
(RW) master functions. Violations of this relationship 
are necessarily due to numerical errors and/or incorrect 
initial conditions. 

We will refer to errors seeded by the initial conditions 
as "junk". One type of junk either propagates off the 
computational domain or decays away. We collectively 
refer to such junk radiation, junk quasi-normal ringing, 
and junk Price tails as dynamical junk. The key ob- 
servation of this paper is that trivial initial conditions 
may also give rise to a static distributional junk solu- 
tion "I 7 jost, which we refer to as Jost junk. In terms of 
the "Schrodinger operator" H = —d% + V, a Jost solu- 
tion satisfies H^f ost = v 2 ^f ost , with ^f ost ~ exp(±wx) 



a, b 


Endpoint of computational domain [a, b]. 


Sl, Sr 


Number of subdomains to left and right of particle. 


N 


Number of points on each subdomain. 


T,S 


Smoothing parameters introduced in Eq. (pjl. 


At,t F 


Timestep and final time. 


M = 1 


Schwarzschild mass parameter. 


m p = 1 


Particle mass. 



TABLE I: Basic set of parameters for a numerical 
simulation. This set is not complete, but in what follows we 
often refer to these variables. For all our simulations M = 
1 = m p , where the choice m p — 1 is equivalent to working 
with per-particle-mass perturbations $/m p . 



as x — > oo [2J]. In this paper, we are exclusively inter- 
ested in "zero-energy" Jost solutions for which v = 0, in 
which case ^jost docs not behave exponentially at infinity 
(see below). Therefore, in what follows a Jost function 
satisfies a "zero-energy" , time- independent, Schrodinger 
equation (— fl£ + F) v Pj os t = to the left and right of the 
particle, and, as it turns out, is discontinuous at the par- 
ticle location. We find that ^Pjost has a non-negligible 
effect in the wave-zone, yet is often small enough to be 
buried into the 0(r _1 ) error associated with a waveform 
"read-off" in the far-field. 

We will adopt trivial initial conditions throughout, but 
allow for modified "smoothed" source terms according to 
the aforementioned description ([2]). Our chief goal is to 
study the properties of the numerical solutions computed 
with and without smoothed source terms, especially in 
the context of the Jost solution. To carry out numeri- 
cal simulations, we have primarily used the nodal Legen- 
dre discontinuous Galerkin method described in Ref. [15| , 
and further details of this method will not be given here. 
In addition, some of our results have either been ob- 
tained or independently verified with a nodal Chebyshev 
method (similar to the one described in Rcfs. [2^, [26|). 
which also features multiple subdomains and upwinding. 
Our nodal Chebyshev method treats the jump discon- 
tinuities at the par ticle location in the same fashion as 
outlined in Ref. [l5| for the nodal dG method. Both our 
dG and Chebyshev methods solve a first order system 
representing (JTJ) . Namely, 

d x V = - II (3a) 
d x ll = /3«c\II - (dx/dQ^deUdx/dt)- 1 *] 

a A $ = a^*)-a f n + j 2 *(f-€ P ) J (3c) 

where the time-space coordinates (A, £) are adapted to 
the particle history (the particle location £ = £ p remains 
fixed in this system). Eq. (|3k) defines II, the variable 1 
<I> = d^, and Ref. [f5[ relates the A-dependent jump 



In our approach, from all fields we explicitly remove delta 
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FIG. 1: Dependence of C on smoothing parameters. 
We have empirically determined that \Cl \ = \ = \Cr\ for an 
impulsive start-up, corresponding to C = 1 at the leftmost 
point. The parameter 8 is different for each r; 8 — 2 for 
r = 10 and 8 = 0.0058 for r = 150. 



terms J1.2 to the sources in ([T]). Most of this paper con- 
siders circular orbits, for which A = t, £ = x, and the shift 
vector /3* = 0. We often refer to the variables II and $ 
below, and for circular orbits these are —dt 1 ^ and d x ^, 
respectively. Throughout the paper, we make reference 
to the parameters listed in Table |U 

This paper is organized as follows. Section UU focuses 
on the Jost solution, from both empirical and analyti- 
cal standpoints. Here we present analytic formulas for 
Jost solutions and compare them with numerical results. 
Section lnTl considcrs several practical consequences of im- 
pulsive start-up for EMRB modeling with circular orbits: 
violation of the axial consistency condition, contamina- 
tion of waveform luminosities, and influence on self-force 
measurements. This section also gives a preliminary re- 
port on consequences for eccentric orbits. Concluding re- 
marks are given in Sec. IIV| where we touch upon finite- 
difference methods while discussing the universality of 
our results. Longer calculations appear in the Appendix. 



II. JOST SOLUTION 

To better explain the origin of the Jost junk solution, 
we first consider a toy model: the ordinary 1+1 wave 
equation with distributional forcing. We then examine 
the Jost solution for the master wave equations, with a 
forcing determined by a circular orbit. 



function terms arising from the distributional inhomogeneity. 
Therefore, <E> = does not hold in the sense of distribu- 

tions. More precisely, in the case of circular orbits, our <E> is 
9**- [[<K]]8(x - x p ). 



A. Forced 1+1 wave equation 

For a fixed velocity v obeying \v\ < 1, we consider the 
model 

-d 2 ^ + flg* = G(t)8(x - vt) + F(t)6'(x -vt) 
G(t) = cost = -iF(t). 

Ref. [HI has shown that 

\&(t, x) = —§ sinz? + \\^ 2 [v + sgn(x - vt)] costf 

(5) 

$ = 7 (t — xv — \x — vt\) 

is an exact particular solution to Here 7 = (1 — 

v 2 )- 1 / 2 is the usual relativistic factor. For this model, 
junk radiation propagates off the computational domain 
with speeds ±1. However, when numerically solving this 
equation subject to (incorrect) trivial initial conditions, 
we observe that the numerical solution no longer con- 
verges to the particular solution. For simulations involv- 
ing ([3]), we have used the dG method with (cf. Table [I]) 
a = -100, b = 100, S L = 10, S R = 10, N = 27, and 
At = 0.01. To compute errors relative to the exact solu- 
tion, we have first interpolated onto a uniformly spaced 
x-grid with 5121 points. Furthermore, to better model 
the circular orbit scenario for EMRBs, we have taken 
v = 0. 

With the exact solution used to generate initial con- 
ditions at t = 0, the nodal dG method exhibits spectral 
convergence throughout the computational domain (and 
for all fields with the wave equation treated as a first or- 
der system) (l5j . However, with trivial initial conditions, 
only the corresponding numerical derivatives, n numer j C ai 
and (I 1 numerical, converge to the correct values, whereas 
^numerical itself is off by a constant value on each subdo- 
main. Let us write 

^numerical = (#L + C L ) 0(-x) + + C R ) Q(x), (6) 

where &(x) is the Heaviside function and the exact solu- 
tion from ([5]) is 

^ l = — \ sin(i + x) — |i cos(t + x) 
* r = — g sin(i — x) + |-i cos(t — x) . 

We introduce the time-independent 1+1 Jost junk solu- 
tion 

*jost = C L Q{-x) + C R Q{x), (8) 
in order to express the numerical solution as ^numerical = 

^exact + *Jost- 

We examine the dependence of C = |Cx| + \Cr\ on the 
smoothing parameters (r, <5), defined analogously to those 
in ([2]) , but here introduced to smooth our toy source term 
cost8{x) + \ cost5'{x). We restrict the parameter space 
by first choosing r, and then finding the smallest <5 such 
that £[erf(\A(i - t Q - r/2) + 1] is less than 10 -16 when 
t = and greater than 1 — 10~ 16 when t — r. These 
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requirements ensure that the start-up phase is smooth to 
machine precision, while providing the most gradual rate 
at which the distributional source terms are turned-on. 
Figure [1] shows that the troublesome constant term is 
arbitrarily well suppressed by the smoothing procedure. 
However, we find that the value of C remains fixed when 
varying the timestep. The final run time for each data 
point in the plot is tp = r + 150. No essential difference 
exists between the v = and u ^ cases, except that 
for the latter case we must ensure that the particle does 
not get too close to the boundary. Let ^smooth represent 
'J' numerical obtained with smoothing, and ^impulsive rep- 
resent ^numerical obtained without smoothing. Then we 
have shown <£ smooth — ^ exact, so that 



Jost 



DOth 



(9) 



is another expression for the Jost solution, valid up to 
method error. In the next subsection we consider this 
expression in the context of master wave equations. 



B. Master wave equations 

The first numerical experiment in this subsection in- 
volves the axial sector with 



V 



axial 



(r) 



f(r) 



r 



(10) 



in (TTI). and assumes CPM source terms (see the appendix 
of [l5| for the precise expressions). To empirically ver- 
ify that an impulsive start-up also leads to a Jost solu- 
tion in this setting, we will form and plot the expres- 
sion (J9)), using the Chebyshev method. Later on, we will 
give analytic expressions for static Jost solutions. The 
experiment enforces Sommcrfcld boundary conditions at 
the left physical boundary, and radiation outer boundary 
conditions [l5l l27j on the right boundary. Our smoothing 
parameters are r = 150 and S = 0.0058. We compute the 
(£, to) = (3,2) metric perturbations for a particle in cir- 
cular orbit initially at (r, <f) = (7.9456,0). Other param- 
eters (cf. Table HJ) are a ~ -202.16, b = 60 + 21og(29) ~ 
66.73, S L = 30, S R = 8, N = 26, At ~ 0.03, and 
tp = 600. Figure [5] shows the result. The plots sug- 
gest that the Jost junk solution affects VP 
spatial derivatives. 

For both axial and polar perturbations generated by 
circular orbits, we now present the analytic form of the 
Jost solution, suppressing throughout the analysis both 
orbital i and azimuthal m indices. For circular orbits we 
have observed empirically that the Jost junk solution can 
be written as 



CPM anH its 
impulsive ana ltS 



axial/polar 
Jost 



axial/polar ^ / 

C L v L lv 9(- 



\ . s~i ax 
X) + Crv r 



il/polar 



G(x), 



(11) 



^impulsive ~ *s 
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FIG. 2: Difference between smoothly and impulsively 
started CPM fields. Here I = 3, m = 2, and the snapshot 
is taken at t = 600. 



where Cl and Cr are complex constants. The functions 
v l?r P ° laI satisfy a Schrodinger equation Hv = defined 



by the operator 



^axial/polar q2 _|_ ^axial/pola 



(12) 



where y axial is given in Eq. (TTJJ)) and, in terms of n = 
\{l -!){(> + 2), 

2 /(0 „ 



V rpolar (r) 



(nr + 3M) 5 



n 1 + n + 



3M\ 9M 2 



M 



(13) 

The functions l , axial /p°i ar ga^fy ^ ne Schrodinger equation 
to the left of the particle, and the functions w ^ xlal / Polal 
the equation to the right. The relevant solutions to Hv = 
decay either as r — > 2M + or r — > oo. 

We derive expressions for all four functions w a ^ al /P olal 
in the Appendix, adopting the dimensionless radius p = 
(2M)~ 1 r as the basic variable. Here we record the set of 
axial functions, 

vl™\p) = p-^F x (l + j + l,t - j + 1; 1; (p - I)/ p) 

(14a) 

^ xial (p) - p-'tF^l + J +l,t- 3 +l-2{l+ 1); p" 1 ), 

(14b) 

where for gravitational perturbations the spin j = 2. Ev- 
idently, up to transformations of the dependent and inde- 



pendent variables, the equation H 



axial , 



is the hyper- 



geometric equation. The equation H polm 'v = involves 
an extra regular singular point, and its normal form is 
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a particular realization of the Hcvm equation. Neverthe- 
less, by exploiting certain intertwining relations between 
the polar and axial master functions (28[ , we are likewise 
able to express v^r in terms of the classical Gauss- 
hypcrgcomctric function 2-Fi- The Appendix gives fur- 
ther details. 

To complete our analytic expressions for the Jost solu- 
tions, we still must determine Cl and Cr. Our notation 
for a time-dependent jump is, for example, 

[[*]] (t) = lim [*(t, r p (t) + e) - *(t, r p (t) - e)] 

= lim [*(t,rp + e)-*(t,r p -c)], (15) 

with the last equality holding for a circular orbit. As 
derived in Ref. [l5j], for a circular orbit the analytical 
jump determined by Eq. ([!]) is 



[[^analytic]] (t) 



(16) 



where the subscript "p" indicates evaluation at the par- 
ticle location. For trivial initial data (that is 3/ = 0) 
this jump will in general not be satisfied at t = 0. We 
find empirically that the jump in 'J'jost exactly cancels 
[[* analytic]] (0), while the jump in c^Jost is zero. The 
system of equations used to determine our constants is 
therefore 



vn(r P )C R - v L (r p )C L = - 
v' R {r P )C R - v' L (r p )C L = 



F(0,r p ) 
fv 



(17) 



which has solution 



Cr 



F(0,r p ) 



fv 



vrv l - VLV R 



Cl = Cr 



(18) 



Recall that ^jost may be numerically approximated as 
^impulsive - * smooth [cf. Eq. ©] . Figure [3] depicts the 
relative error |(* Jost - (^impulsive - *smooth))/*Jost | for 
£ = 3 perturbations, with ^jost given by (fTTj) . To gen- 
erate this figure, we have used nearly the same set-up 
as described for Fig. [21 but with the outer boundary 
b = 240 + 21og(119) and final time t F = 3100. 



C. Jost solution and radiation boundary conditions 

We wish to examine the extent to which the right ana- 
lytic Jost solutions w ^ xlal /P° lal satisfy radiation boundary 
conditions based on Laplace convolution [13, H^, as these 
are boundary conditions adopted for our numerical sim- 
ulations. Unfortunately, for blackhole perturbations the 
issue would seem difficult to address analytically. There- 
fore, we consider the analogous issue for the flatspace 
radial wave equation. 



Consider a flatspace multipole solution 
r _1 \I/(i, r)Yg m (9, (j>) to the ordinary 3+1 wave equation, 
and assume the multipole is initially of compact support 
in radius r. Exact non-reflecting boundary conditions 
relative to a sufficiently large outer boundary radius b 
then take the form ^9f 





9*\ 


\~dt ' 





r—b 



b 2 ^ 



exp (b^hjit - *'))*(*', b)dt'. 



(19) 



Here {kej : j = 1, ...,£} are the roots of the modi- 
fied cylindrical Bcssel function K( + i/ 2 (x), also known as 
MacDonald's function. All kij lie in the left-half plane. 
Moreover, the scaled roots ke t j/(£ + 1/2) accumulate on 
a fixed transcendental curve as I grows [13, HH, so the 
exponentials exp (b ktjt) tend to decay more quickly 
in time t > for larger I. 

For the flatspace setting at hand, the Jost solution sat- 
isfies 



£(£ + 1) 



0. 



(20) 



and two appropriate linearly independent solutions are 
the following: 



vr(t) 



(21) 



We therefore examine to what extent vr{t) satisfies ([T9 
Straightforward calculation yields 



m 



r—b 



= -6 -1 »k(6) ^2 ex P 



1 1 f* 
-- Y^h,j / exp^- 1 

.7 = 1 J ° 



(22) 



k t j(t-t'))vR{b)dt'. 



The function vr(t) does not satisfy the non- reflecting 
condition (|19p ; however, the violation of (fl9|) decays ex- 
ponentially fast. For blackhole perturbations we likewise 
expect that i;^ polar (/9) violates our radiation bound- 
ary conditions only by exponentially decaying terms, and 
have seen some evidence of this behavior in our numerical 
simulations. 

We have also observed persistent junk solutions when 
adopting the Sommerfeld condition at the outer bound- 
ary b along with impulsive start-up. We differentiate be- 
tween two scenarios: the first involving a detector which 
is not in causal contact with the outer boundary b dur- 
ing the simulation, and a second with the detector lo- 
cated at b. For the first scenario, the static junk solution 
which develops and persists around the detector is pre- 
cisely ^jost- For the second, we also observe a persistent 
junk solution, but one which is distorted from \tjost in a 
boundary layer near b. Such distortion presumably arises 
since ^j st satisfies the outer Sommerfeld condition only 
up to an 0(r - ^ -1 ) error term. 
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FIG. 3: Comparison between analytic and numerical Jost solutions. CPM and ZM modes respectively correspond to 
{£, m) = (3, 2) and (£, m) = (3, 1), (3, 3). 



III. CONSEQUENCES OF IMPULSIVE 
STARTING CONDITIONS 

A. Inconsistent modeling of the axial sector 

Axial perturbations are described by either the 
Cunningham-Pricc-Moncrief master function iJ/ CPM or 
the Reggc- Wheeler master function v]/ RW . Both solve the 
generic wave equation (JXJ) with potential (JTUJ). However, 
the wave equations for v[/ CPM and v]/ RW have different dis- 
tributional source terms [13, HE H(| ■ As shown in [3(| , 
these master functions obey 



1> 



nw 



IttCPM 
2 11 



= 0, 



^ r P (t), 



(23) 



and we refer to this formula as the axial consistency 

condition. For circular orbits this condition becomes 
vjjRW _ i^cpm _ Qj r ^ ^ We nQW numer i ca iiy 

examine the extent to which the axial consistency condi- 
tion is violated when the master functions ij/ RW < CPM are 
obtained with and without smoothing. 

For all experiments we again enforce Sommcrfcld 
boundary conditions at the left physical boundary, and 
radiation outer boundary conditions on the right bound- 
ary. Now our smoothing parameters are to = 0, r = 100, 
and S = 0.05. We compute the (£, m) = (2, 1) met- 
ric perturbations for a particle in circular orbit initially 
at (r,<t>) = (7.9456,0). Other parameters (cf. Table Q} 
are a = -200, b = 30 + 21og(14) ~ 35.28, S L = 22, 
S R = 3, N = 31, At = 0.01, and t F = 800. We first 
plot |* RW + in CPM | at various times. The left panels 
in Fig. |4] show results with smoothing. Although the 
consistency condition is initially violated, the expression 
eventually relaxes to a small value once the dynamical 
junk has propagated off the domain. The right panels in 
Fig.|4]show result without smoothing. Even at late times 
violation in the axial consistency condition is now evi- 
dent. The plots in Fig.[5]dcpict |* RW + in CPM | recorded 



as a time series at x = —200. The plot for smooth 
start-up indicates that quasinormal ringing and Price de- 
cay tails characterize the late-stage dynamical junk, al- 
though this ringing is suppressed with more smoothing 
(c. g. with r = 150, 6 = 0.0058). The plot for impul- 
sive start-up suggests that a static Jost junk solution 
*^ui si ve " *25Lh P ersists indefinitely (n CPM should 
be unaffected by a similar Jost solution in \I> CPM ). 



B. Contamination of waveforms 

For a given (£, m) multipole either read off at a finite 
radius or measured at null infinity through an approxi- 
mate extraction, we can apply standard formulas to esti- 
mate the energy and angular momentum carried away by 
the gravitational waves. We continue to work with the 
axial perturbations, with formulas featuring only CPM 
and RW masterfunctions. The luminosity expressions are 
the following: Q3, 



E, 



CPM 



r CPM 



RW 



E, 



f RW _ 



1 (£ + 


2) 


I^CPMh 
| Ivn 


64tt {£ - 


2) 


im (£ + 


2) 


vjyCPM^ 


64tt (£ - 


2) 


1 (£ + 


2) 


I^RWl 2 

| Ivn \ 


16tt (£ - 


2) 


im (£ + 


2) 


/ 


16tt (£ - 


2) 



In the distant wave-zone we expect Ef^ M 

'CPM 



(24a) 



(24b) 



EJ™ and 



L^ m = L R w. However, Sec. HIS has shown that im- 
pulsive start-up can result in violation of the axial consis- 
tency condition (|23p , and such violation in turn results in 
discrepancies between the above luminosity formulas. As 
seen in Sec. Ill Bl whether simulations are based on vj/ CPM 
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FIG. 4: Snapshots of |^ rw + |n CPM | with and without smoothing. The left three panels correspond to smooth start-up 
and the right three to impulsive start-up. The times at the far left correspond to both sets of panels. VJ> RW ; s of order 10~ 2 
near r v . 



or vj/Kvv^ an impulsive start-up generates a Jost junk so- 
lution, even at long distances from the source. Although 
dynamical junk is also present, its effect is negligible in 
the wave-zone at late times. 

Table [XT] collects summed luminosities for (I, m) = 
(2,±1) waveforms. The top set of numbers arc un- 
averaged and recorded at time tp = 2750, while the 
bottom set have been averaged between t = 2500 and 
t F = 2500 + 4T , where T = 2tt P 3 / 2 ~ 140.7246. 
Other parameters (cf. Tabic Q} arc a ~ —190.34, b = 
1000 + 21og(499) ~ 1012.43, S L = 30, S R = 150, 
iV = 26, and At = 0.038. For smoothing we use 
t = 150 and 6 = 0.0058. For circular orbits we ex- 
pect (Osmooth) = Qsmooth, where brackets denote time 
averaging for a generic luminosity Q. Relative errors are 
computed by 

A, | ^smooth ^impulsive) /oc:^ 

terror 7~p~ j • \ ) 

| smooth | 

For the CPM luminosities computed with smoothing, 
time averaging has little effect. However, it does en- 
hance the accuracy of the RW luminosities computed 



with smoothing. Indeed, inspection of the bottom sec- 
tion of Table |TT] shows that the CPM and RW entries in 
the Qsmooth column are in excellent agreement. 

Relative to the true luminosity which would be 
recorded at null infinity, even the exact E CPM read off at 
r = 1000 would have an 0(r _1 ) error, but here we have 
viewed the read-off value as the true one. Because £; CPM 
is unaffected by the Jost junk solution, estimates 
error stemming from both the method (here the Cheby- 
shev scheme) and any residual dynamical junk. The 
other luminosities are affected by the Jost junk solution; 
however, as shown in the Appendix, errors which stem 
from the Jost solution decay faster than 1/r. Therefore, 
these errors should be smaller than the 0(r _1 ) errors as- 
sociated with using the read-off luminosities as approxi- 
mations to the ones at null infinity. 



C. Self- force measurements 

Over long times the influence of the metric pertur- 
bations on the particle orbit should significantly affect 
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FIG. 5: Time series ati = —200 for |>I> rw + In CPM | with and without smoothing. * rw is of order 10 4 at x = -200. 



Q 


Qsmooth 


Qimpulsivo 


Qcrror 


gCPM 
^CPM 

i RW 


8.17530620 x 10~ 7 
8.17530652 x 10~ 7 

1.83102415 x 10~ 5 + i3. 24326408 x 10~ 14 
1.83047467 x 10 -5 - i2. 16502183 x 10~ 8 


8.17530623 X 10~ 7 
8.18248752 x 10~ 7 

1.82972897 X 10~ 5 - il. 28610911 x 10~ s 
1.66825388 x 10~ 5 +i8. 14152318 x 10~ 7 


3.4668 x 10~ 9 
8.7838 X 10~ 4 
9.9685 X 10~ 4 
9.9693 x 10~ 2 


(B OPM ) 
(E KW ) 

(i OPM ) 
<i RW ) 


8.17530620 x 10 -7 
8.17530617 x 10~ 7 

1.83102416 x 10~ 5 - il. 40467882 x 10~ 15 
1.83102415 x 10~ 5 +i4. 13269715 x 10~ 13 


8.17530620 x 10~ 7 
8.17531431 x 10~ 7 

1.83102416 X 10~ 5 +i3. 49294212 x 10~ 14 
1.82927679 X 10~ 5 + i7.05636411 x 10~ 9 


2.8376 X 10~ 10 
9.9661 x 10~ 7 
2.0738 X 10~ 9 
1.0292 x 10~ 3 



TABLE II: I = 2 luminosities recorded at r = 1000. Entries result from addition of m = 1 and m — — 1 luminosities, 
and they correspond to a circular orbit with (r, cj>) = (7.9456, 0) initially. Qorror as been computed with more precision than 
reported for the table entries. 



the gravitational waveform [3l| . and realistic waveform 
computations will therefore need to include this influ- 
ence. Incorporation of self-force effects constitutes one 
approach towards modeling this influence. Because the 
metric perturbations are discontinuous at the particle, 
self-force calculations typically require a regularization 
technique. In the Regge- Wheeler gauge no regulariza- 
tion procedure exists f or g eneric orbits; however, direct 
field-regularization [32l . |33| | seems promising. For the re- 
stricted case of circular orbits, Detweiler has shown how 
to directly calculate certain gauge invariant quantities in 
the RW gauge without regularization [34 ]. Detweiler's 
approach obtains the energy luminosity Eqw and an- 
gular momentum luminosity Low associated with waves 
escaping to null infinity and down the black hole through 



local self-force calculations performed at the particle, 



E --—u a u^ — 
p ~ 2u* 



at 



L, 



1 

2u* 



dh 



a/3 



(26) 



where the perturbation h a p is reconstructed from <F 
and its derivatives [3f|. Equations (|26|) hold for each 
(£, m) mode of the metric perturbation. For pertur- 
bations described by the CPM masterfuntion and with 
the Regge- Wheeler gauge, the non-zero contributions (for 
each mode) involve 



dh t 



dt 
dh t4 



9 2 * 
dtdr 

a* 

dr 



<9<F 



(27) 



in a source free region. Here X^ and X^ are axial vector 
and tensor spherical harmonics [3fj| . When numerically 
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forming these expressions, we replace dt*& and d r ^ by — II 
and Only when evaluated at the particle location 

will E p and L p be related to Eqw and Lqw- 

We now fix r = 100 and S = 0.014 to achieve a 
smooth start-up, run to the final time tp = 800, and pick 
At = 0.005. Other parameters are the same as those in 
Sec. IIIIAI We compute E p and L p for (£,m) = (2,±1) 
perturbations. Because E p is computed with time deriva- 
tives of \I/ CPM ; the static Jost junk solution does not im- 
pact its measurement. We therefore expect that 

^p(*impulsive) — ^(^smooth)- ( 28 ) 

However, an impulsive start-up appears to generate more 
dynamical junk at late times. Figure [5] depicts E p , 
recorded as a time series, for both impulsive and smooth 
start-ups. A separate experiment based on waveform 
read-off near the blackhole and waveform extraction at 
the outer boundary determines that the energy carried 
away by the gravitational waves is Eqw — 8.3163 x 10~ 7 . 
The relative errors in the left panel of Fig. [5] arc com- 
puted as \E p — Eqwj\/Egwj and are limited by the ac- 
curacy of Eqw- We therefore do not expect agreement 
beyond a relative error of 10~ 5 , although clearly such er- 
ror will settle to a constant value. The time series for 
both the impulsive and smooth start-up exhibit large os- 
cillations which persist until about t = 400. However, 
beyond t = 400 the impulsive start-up series shows larger 
oscillations. 

L p depends on both \[/ CPM and its spatial derivative 
$ CPM , whence the Jost junk solution will impact its 
self force measurement. With smoothing, the time series 
plot for Lp looks similar to one for E p in Fig. [6j and is 
not shown. We note that our self-force L p measurement 
agrees with a separate experiment which finds that the 
angular momentum carried away by gravitational waves 
is L G w — 1.8626 x 10~ 5 . Figure [7] shows that L p is 
typically discontinuous at the particle for an impulsive 
start-up. Even with an impulsive start-up, the L p mea- 
surement yields the correct value when averaged over an 
orbital period X^, and it is continuous across the particle 
(with the correct value) when the particle returns to its 
initial orbital angle. 

These phenomena are a consequence of the axial Jost 
junk solution (fTT]) . For t fixed, Eq. (|2"6")l shows that 
L p (^) depends linearly on "P. Therefore, L p (tyj™ st + 

Smooth) = £p(*j£t) + 4 (Smooth), SO WC Can foCUS 

on ip^j^ijt) alone. The expressions ([15} for Cl,r are 
linear in F(0,r p ), which is in turn proportional to the 
conjugate of an axial vector spherical harmonic X^ [30l | . 
Motivated by this observation, we "factor off" the con- 
jugate, writing *^™ t = 7j e (x)X^ m ((j) ), where O is the 
particle's initial orbital angle and r\i (x) is a real discon- 
tinuous function solely of x. The expression (|26[) for L p 
involves dht^/dcj), which by (|27p is proportional to X^. 
In the equatorial plane X^™ = d^X^ m = imXj™, and we 

conclude that L p (¥$£t) = im^{x)X^{<j> )X^ v {t)), 
where &(x) is a real discontinuous function solely of x. 



Therefore, when the particle returns to its initial position 
(that is, when 4> p (t) = (f>o), the value of ip(^j^ t ) is pure 

imaginary and ip(*j™ t ) + ip(*j st") = 0- For pertur- 
bations generated by a particle in circular orbit, we have 
seen that *f™ pulsive s + *f™ ooth to high accuracy. 

Combination of this expression and the above arguments 
for axial perturbations then gives 

E 4(*impul S ive) * £ L p (¥™ ooth ), (29) 

\m\<e \m\<£ 

when <j) p (t) = 4> - Moreover, one finds (i p j™ t ) ) = 
for time averaging over an orbital period T^. 



D. Consequences for eccentric orbits: preliminary 
results 

This section considers a particle in the eccentric or- 
bit described in Section IV. B. 2 of [l5j]. In the notations 
of that reference the orbit's eccentricity and semi-latus 
rectum are (e = 0.76412402, p = 8.75456059), and we 
choose x = 0.2 and <f> = ir/4 to fix the particle's ini- 
tial position. We simulate the resulting (£, m) = (2, 1) 
perturbation with (cf. Tabic H} a = -200, b = 1012.43, 
S L = 22, Sr = 100, N = 31, At = 0.02, and t F = 3000. 
We again take r = 150, S — 0.0058 as the smoothing 
parameters. Since e ^ 0, we use a coordinate transfor- 
mation to keep the particle at a fixed location between 
subdomains (see [l5[ for details) . Before making compar- 
isons, we first interpolate all fields onto a uniform x-grid 
(tortoise coordinate) with 6063 points. 

Fig. [8] shows the difference between fields for smooth 
and impulsive start-ups. The two numerical solutions are 
clearly different, although for the case of eccentric or- 
bits we have no analytical understanding of the resulting 
"junk solution" 2 presumably seeded by impulsive start- 
up. Empirically, we find that this solution satisfies 

[[*ju„k]](t) = -[[*a„alytic]](0) (30a) 
P>ju„k]]W=0 (30b) 
pjunk]] (t) = 0, (30c) 

where [[* analytic]] (*) = f P (t)F(tr p (t)) / (f*(t) - fj»(t)) in 
terms of f p (t) = f(r p (t)). See [l5|| for a derivation of the 
analytical jump. These time independent jump condi- 
tions are the same as for the circular orbit ^jost solution. 
With our choice of numerical parameters the axial con- 
sistency condition is satisfied to better than a 1 x 10 -6 
relative error throughout the entire domain for a smooth 



2 At present, we are uncertain if the generated junk solution fulfills 
the formal definition of a Jost solution. Thus, in the context of 
eccentric orbits we simply refer to the persistent solution as the 
"junk solution". 
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FIG. 6: E p time series FOR summation OF £ — 2 AND m = ±1 MODES. In the right panel the curve corresponding to 
impulsive start-up has the larger amplitude (due to small fluctuations this curve does not appear dashed as indicated in the 
legend) . 
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FIG. 7: t = 800 SNAPSHOT OF REAL part OF L p FOR 1 = 2 AND m = 1. The particle is located at the interface between the 
two subdomains. 



start-up. For an impulsive start-up the condition is vio- 
lated to the order of the solution itself. We conclude that, 
as for circular orbits, the junk solution generated by an 
impulsive start-up leads to an inconsistent modeling of 
the axial sector. 

Table IIIII collects energy and angular momentum lu- 
minosities. These luminosities have been averaged from 
t = 1700 to t F = 1700 + 4T r , where T r ~ 780.6256 is 
the radial period (see [HI for further details). Unlike 
the circular orbit case, the discrepancy between wave- 
forms corresponding to smoothly and impulsively started 
fields may be larger than usual 0(1/ r) error associated 



with read-off at a finite radial location rather than in- 
finity. Moreover, the junk solution would seem deter- 
mined by the initial orbital parameters. Indeed, the val- 
ues Qimpuisivo and errors quoted in our table strongly 
depend upon such choices. 



IV. CONCLUSIONS 

A number of time-domain methods exist for solv- 
ing Eq. ([1]) as an initial boundary value pro blem, in- 
cluding those described in @, Q3, [U H [H, [H, 
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FIG. 8: Difference between CPM fields with and without smoothing for an eccentric orbit. Here we plot both 
real (dashed) and imaginary (solid) parts at ip = 3000. 



Q 


Qsmooth 


Qimpulsivc 


Qerror 


(B 2 c r M > 


1.559917 x 10" 4 
1.153983 x 10~ 6 
1.153983 x 10~ 6 


1.559484 x 10~ 4 
1.236758 x 10~ 6 
1.872073 x 10~ 6 


2.775789 x 10~ 4 
7.172983 x 10~ 2 
6.222709 x 10" 1 


<B 2 c r M ) + 


1.571457 x 10~ 4 


1.571852 x 10" 4 


2.512000 x 10~ 4 


Rc<i°™) 
R.c(L™) 


2.078556 x 10~ 3 
1.441737 x 10~ 5 
1.441749 x 10~ 5 


2.076811 x 10~ 3 
1.537876 x 10~ 5 
1.662726 x 10~ 5 


8.395251 x 10~ 4 
6.668276 x 10~ 2 
1.532701 X 10 _1 


Rc<i°™> + Rc(L™> 


2.092973 x 10~ 3 


2.092190 x 10~ 3 


3.744004 x 10~ 4 



TABLE III: £ — 2 luminosities for a particle with an orbit given by (e = 0.76412402, p = 8.75456059). Entries result 
from the addition of \m\ and — |m| luminosities. 



These methods vary in both the underlying numerical 
scheme (e.g. finite difference, finite element, pseudospec- 
tral, and spectral) as well as their treatment of the dis- 
tributional source terms (e.g. Gaussian representation, 
analytic integration, domain matching). Numerical sim- 
ulation of metric perturbations may also involve other 
choices (e.g. gauge, number of spatial dimensions, choice 
of numerical variables). Moreover, similar time-domain 
methods exist for solving the forced Teukolsky equation 
describing particle-driven perturbations of the Kerr ge- 



ometry (see for example Refs. (38H40j ). For all of these 
methods, the issue of impulsive start-up would seem per- 
tinent, although clearly we cannot examine each method. 
Nevertheless, we now attempt to provide at least partial 
insight into the ubiquity of static junk solutions. 

As mentioned earlier, the results and observations of 
this paper have been independently confirmed with each 
of our two numerical methods: the nodal Legendre dG 
and Chebyshev schemes. However, as these schemes are 
rather similar, we now briefly consider a finite-difference 
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FIG. 9: Difference between smoothly and impulsively 

STARTED FIELDS USING A FINITE-DIFFERENCE METHOD. As in 

Subsection HlBl we consider * CPM for I = 3 and m = 2. The 
bottom plot depicts the relative error between the numerical 
and analytical Jost solutions. 



scheme for solving ([3]), based on fourth, sixth, and eighth 
order stencils for the spatial derivatives. To stabilize 
sixth and eighth order stencils, we have followed Rcf. 



41 



Furthermore, we replace the Dirac delta functions in ([3]) 
by narrow Gaussians. Precisely, for cr = 0.1 wc make the 
replacement 



J(x, t)S(x — x p ) — > J(x, t) 



2na 



exp 



2a 2 



(31) 

for both the J\ and J 2 terms in ([3]). Analytic expressions 
for Ji and J2 are readily computed with Eq. (28) from 
Ref. [15| ■ With essentially the same experimental set-up 
described in Subsection III Bt we repeat that experiment 
using 4000 points and sixth order spatial differences. The 
results, shown in Fig. [SI clearly indicate the presence of a 
static "Jost junk" solution. A larger choice for a gives rise 
to a rounder transition near the particle. However, the 
following shows that such contamination is not a generic 
feature. For circular orbits, our system ([3]) becomes 



a t tf = -n 

d t U = -d x <f> + V(r)* + Ji6(x 
d t $ = -djl + J 2 5(x - x p ), 



(32) 



so that the system formally becomes 
d t $ = -n 

d t U = -d x $ + V(r)% + JiS(x - x p ) + J 3 5'(x - x p ) 

dt$ = -aji, 

(34) 

where J3 = [[\&]] = F(t, r p )/ f p . If we now replace the 5,5' 
terms in the new system by appropriate Gaussians, then 
we do not observe a persistent Jost junk solution when 
trivial initial conditions are supplied (neither in finite- 
difference nor dG simulations). 

Persistent junk solutions arise from the combination 
of inconsistent initial data and the distributional forcing 
terms which define the EMRB model. In particular, we 
observe that development of a Jost junk solution depends 
on how the distributional forcing is treated rather than 
the underlying numerical method. Therefore, whether or 
not they contaminate simulations should be considered 
on a case-by-case basis. Domain matching approaches 
which enforce jump conditions without approximation 
(considered in this paper) exhibit a Jost junk solution 
in the absence of smooth start-up. With first order vari- 
ables such approaches correspond to system (|32|) rather 
than (|3~4")) . Treatment of system (f3~4"]l with Gaussian rep- 
resentations for 5,5' exhibits no persistent junk solution, 
although such an approach necessarily introduces large 
method error relative to the exact distributional model 
and features variables with 5-like behavior near the "par- 
ticle" (Gaussian peak). The issue of a static junk solution 
for schemes which discrctize the second order equation 
([T]) deserves further consideration, although, if present, 
then the particular Jost junk solution observed in this 
paper would likely be of relevance. 3 

We have shown that impulsive starting conditions are 
inadequate for time-domain modeling of extreme mass 
ratio binaries. Such conditions result in more dynamical 
junk, evident in self-force calculations, and potentially 
a static Jost junk solution which persists indefinitely. 
Although each effect is small compared to the physi- 
cal solution, such systematic errors will corrupt studies 
which require high accuracy. For example, computation 
of waveforms accurate to second order in the mass ratio 
requires reconstruction of the first order perturbations. 
Since these first order terms act as sources for the wave 
equations describing the second order mastcrfunctions, 
the presence of a Jost junk solution will affect second or- 
der waveforms. Circular orbits far from the massive cen- 
tral object (of potential relevance for the quasi-circular 
phase of inspiral) are similarly impacted by the Jost junk 



where the time-dependent jump factors are J\ — [[^ x ]] 
and J2 = — [[^/]]- Wc introduce a variable $ obeying 

$ = $ - p]]<Kx - x p ), (33) 



For a static solution to have gone unnoticed, it would seem rea- 
sonable to expect decay as either r — ► 2M+ or r — > 00. Such 
solutions will necessarily be discontinuous, and presumably such 
discontinuities could only "hide" at the particle, requirements 
that fix the form of the static solution up to the constants Cl 
and Cr introduced in Section III Bl 
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solution. Eq. (jT8|) indicates that the magnitude of a 
polar-mode static junk solutions does not decay as r p 
becomes large (compare with Eqs. (C5a) and (C6e) from 
[l5l| ) . However, such decay is present in the axial case 
(compare with (C8a) and (C9c) of [Hj|). When studying 
eccentric orbits, errors arising from the persistent junk 
solution appear to corrupt studies requiring even mod- 
est accuracy. Minimization of dynamical and Jost junk 
by smoothing the source terms during start-up will im- 
prove waveform templates and self-force techniques with 
minimal computational and human effort. 
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Appendix: Time-independent master equations 

1. Regge- Wheeler equation 

Subject to the Ansatz that the solution v is time- 
independent and in terms of the dimensionless variable 
p =j2M)~ 1 r, the homogeneous Regge- Wheeler equation 
is [H 



1 



1 



1 



e(£ + i) 



v = 0, (35) 



where k = 1 — j 2 in terms of the spin 3 = 0, 1, 2. For 
gravitational perturbations 3 = 2, but we leave 3 unspec- 
ified for the time being. Expressing the equation in the 
form 

v" + P{p)v' + Q(p)v = 
1 



P{p) = 



Q(p) 



£{£+l)p + K (36) 

pHp-i) ' 

we find that it has regular singular points at 0, 1, and 
00, as well as the associated Riemann-Papperitz symbol 



1 00 

P{ 1+3 -(£ + 1) ; P 
1-3 £ 



> . 



(37) 



-.i+j, 



To obtain the standard normal form, we let v = p" 1 J u, 
so that 



1 00 
Pi -£ + j ;p 
-23 £ + j + l 



(38) 



where u satisfies the hypcrgcometric equation 

p(l - p)u" + [c - (a + b + l)p]v! - abu = 0, (39) 

with a = — £ + 3, b = £ + 3 + 1, and c = 1 + 2j. As one of 
the two linearly independent solutions based at p = 00 
(chosen to be the second), we may take 

MP) =p- e - J - 1 2F 1 (£+ J +l,£-3+l; 2(m); p' 1 )- (40) 

Expressed in terms of the original dependent variable, 
i>2 = p 1+3 U2, this solution is our axial/right solution 



.axial 

J R 



(p) given in (fl4b ). To obtain series so- 



lutions based at 1 which are nevertheless valid on (1, 00), 
we follow Leaver [43| and consider the transformation 
i] = (p — 1) I p. Then with w(r\) = v (1/(1 — f?)), we get 

w" + V{rj)w' + Q(rj)w = 

V{,) = ^L, q^^ 1 ^ 1 ^) (41) 



77(1 - 77) ' 

which has the P-symbol 



r)(l - t/) 2 



1 00 

w = p{ -(^ + 1) 1+3 ;v 

£ 1-3 



(42) 



Now let w = (r/ — l) l y so that 



1 00 
pI l + £ + j ;v 
-(21 + 1) 1 + 1-3 



(43) 



solves 



77(1 - i])y" + [C - {A + B + l)rj\y' - ABy = 0, (44) 



with A 



3 + 1, B = £ + 3+1, and C = 1. Therefore, 



we choose v\{p) = v^ lal (p) given in (fT4k ) as both a first 
linearly independent solution and the axial/left one of 
interest. 



2. Zerilli equation 

In dimensionless form, the time-independent Zerilli 
equation is 



1 



1 



1 

PJ f 
8n 2 (n + l)p 3 + 12n 2 p 2 + 18np + 9 
p 3 (2np + 3) 2 



(45) 



v = 0, 



again where n = \{£ — 1){£ + 2). In standard form, the 
equation is 



v" + P{p)v' + Q{p)v = 0, P(p) 



1 



P(P-1) 



Q(p) = - 



8n 2 (n + l)p 3 + 12n 2 p 2 + 18np + 9 



p 2 (p-l)(2np + 3) 2 



(46) 
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This equation has regular singular points at 0, 1, oo, 
and — 3/(2n), with the following associated pairs of in- 
dicial exponents: {1,1}, {0,0}, {£,-(£ + 1)}, {2,-1}. 
The general second order homogeneous ODE with reg- 
ular singular points at Zq, z\, z 2 , and oo has the form 
y" + R{z)y' + S(z)y = 0, with 



R(z) 
S(z) 



An 



Ap 

z - z 
Bp 

(z-zq) 2 (z-zi) 2 (z-z 2 ) 2 
Co C\ C 2 



Z — Z\ Z — Z 2 

B\ B 2 



(47) 



Z — Z2 



where the Ai, Bi, and Cj are all constants subject to 
Co + C\ + C 2 = and the requirement that for each 
i = 0, 1, 2 at least one member of the triple Ai, Bi, and 
Ci must be nonzero (for otherwise Zi would be a ordinary 
point). By expressing all constants Ai, Bi, Ci except 
Co in terms of the indicial exponents {{A^A^} : k — 
0, 1, 2, oo}, we find 



R(z) = 
S(z) = 



1 — Ao — Ao 1 — Ax — Xi 1 — A2 — X- 2 



z- z 
AqAq 



Z — Z\ 

AiAi 



2 — 22 



\2^2 



(z - z ) 2 (z-zi) 2 (z~z 2 ) 2 
AooA'oo — XqX'q— A1A1— A2A2 
(z - zi)(z - z 2 ) 
C (z - zi)(z - z 2 ) 



(48) 



(z - z )(z - Zx){z - z 2 ) ' 

Here —Co is the accessory parameter (44| . and the gener- 
alized Riemann scheme [44| associated with the equation 
is 



1111 

z zi z 2 00 ;z 
Ao Ai A2 Aoo ; —Co 
Ao Ai X' 2 A'oo 



(49) 



The notation is similar to the P-symbol, but also indi- 
cates the type of singular points in the first row (regular 
singular points are indicated by a 1). Wc find the scheme 



111 1 

1 -3/(2n) 00 ;p 

1 2 -(^+1) ;0 
10-1 £ 



(50) 



for the specific case of the timc-indcpcndent Zcrilli equa- 
tion ([4"5]) . 

Upon transforming the ODE specified by (|48|) to nor- 
mal form, we find the new accessory parameter 



9 = -Co 



A (Ai-l) + Ai(Ao-l) 

z - Zl 
Ao(A 2 -l) + A 2 (Ao-l) 

z — z 2 



(51) 



as well as the transformed scheme 



1 


1 


1 




1 






zo 




Z2 




00 




;z 











Aoo " 


1- A + 


AH 


-X 2 ;q 


-A 


Ai-Ai 


x 2 — x 2 


A'oo 


f A + 


Ai- 


f A 2 



(52) 

With the assumptions zq = and z± = 1, this scheme 
corresponds to the Heun equation G"+P(z)G'+Q(z)G = 
in normal form, where 



P(z) = - 



Q(z) 



l+a+b-c-d 



z-1 
(ib 



z — z 2 



qz 2 



(53) 



(z-l)(z-z 2 ) z(z- l)(z- z 2 )' 
Here the transformed schciuo 



1 

zi 




1 

00 : z 



a ;q 



1 — c l — o? c + d — a — b b 



(54) 



is expressed in terms of the constants a, b, c, and d which 
may be related to the above exponent pairs {A&, A' fc }. The 
normal form of (|4"5")) then has the scheme 



111 1 

1 -3/(2n) 00 

2-1 

-3 ^ + 3 



P 

1 - 4n/3 



(55) 



While the preceding analysis both addresses the struc- 
ture of the time-independent Zcrilli equation and reveals 
the asymptotic behavior of the solutions near any given 
singular point, it does not provide concrete analytical ex- 
pressions for the solutions v£°]^ r considered in the main 
text. To obtain such expressions, we use the intertwining 
operators [28| 



-J) 5* 



2 , ,s 3p-l 
-n(n + 1 + -^-r 11 

3 V ; p 2 {Z + 2np) 



(56) 



Using our earlier solutions v'f a ^(p) to the time- 
independent Regge- Wheeler equation, we then get cor- 
responding solutions Vj j °^ T (p) = D + vf^$(p) to |45|) by 
direct application of D+ and the identity 

ab 



d 
dz 



2 Fi(a,b; c;z) 



2 Fi(a + l,6+l;c+l;z). (57) 



Therefore, we have also expressed the relevant polar so- 
lutions in terms of the Gauss-hypergeometric function 
2 F\. The analysis above then shows that we are likewise 
able to express solutions to a particular instance of the 
Heun equation in terms of hypcrgeometric functions. 
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